############## 
############## 
############## Models
remove(list = ls())

base::library(texreg)
base::library(conflicted)
base::library(tidyverse)
conflict_prefer("filter","dplyr")
base::library(ggplot2)
base::library(dplyr)
base::library(here)
conflict_prefer("here", "here")
base::library(systemfit)
#base::library(sf)
#base::library(sp)
#base::library(spdep)
#base::library(spatialreg)
#base::library(rgdal)
#base::library(utils)
#base::library(geodist)
base::library(stargazer)
#base::library(texreg)
base::library(DataCombine)
#base::library(reshape2)

load(here("Data", "Model_Data.RData"))
glimpse(final_data)

summary(final_data$death_rate)
table(final_data$year)

final_data <- final_data %>% 
  mutate(time_trend = ifelse(year == 2004, 1, NA),
         time_trend = ifelse(year == 2008, 2, time_trend),
         time_trend = ifelse(year == 2012, 3, time_trend),
         time_trend = ifelse(year == 2016, 4, time_trend),
         time_trend = ifelse(year == 2020, 5, time_trend)) 

glimpse(final_data)

eq1_d <- dem_rep ~ death_rate + unemp_rate + 
  log(defl_pcincome) + 
  share_black + share_hisp + share_asian + share_young + share_elder + 
  log1p(share_college) + log1p(share_manuf) +
  log(tot_pop) + as.factor(rural_urban) + state_partiesdiff +
  state + rep_inc + incumbent_cand + lag_dem_rep + time_trend
eq2_d <- abs_rep ~ death_rate + unemp_rate + 
  log(defl_pcincome) + 
  share_black + share_hisp + share_asian + share_young + share_elder + 
  log1p(share_college) + log1p(share_manuf) +
  log(tot_pop) + as.factor(rural_urban) + state_partiesdiff +
  state + rep_inc + incumbent_cand + lag_abs_rep + time_trend
Syst_death <- list(eq1_d, eq2_d)
model_d <- systemfit(Syst_death, method = "SUR", data = final_data)
summary(model_d)

(estia <- cbind(Estimate = coef(model_d), confint(model_d, level = 0.90)))
length(estia[,1])
estia[2,]
eq1ia <- estia[2, ] #death rates - Democrats
estia[75,]
eq2ia <- estia[75, ] #death rates - Abstention

equations <- c("Democratic/Republican", "Abstention/Republican")
results <- rbind(eq1ia, eq2ia)
results

results <- data.frame(cbind(results, equations))
glimpse(results)

results <- results %>% 
  mutate(equations = factor(equations, 
                            levels = c("Democratic/Republican", 
                                       "Abstention/Republican"))) %>% 
  mutate(Estimate = as.numeric(Estimate)) %>% 
  mutate(X5.. = as.numeric(X5..)) %>% 
  mutate(X95.. = as.numeric(X95..))

fig_sur <- results %>% 
  ggplot(aes()) +
  geom_hline(yintercept = 0, colour = "gray37", lty = 2, linewidth = .75) +
  geom_pointrange(aes(x = equations, y = Estimate, ymin = X5..,
                      ymax = X95..), shape = 15, 
                  lwd = 0.7, position = position_dodge(width = .5)) + 
  theme_bw() + 
  ylab("Log-Ratio Form") + xlab(NULL) + #ylim(-0.1, 0.05) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 13),
        axis.text.y = element_text(size = 10),
        axis.title.y = element_text(size = 13),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 14),
        legend.position = "bottom",
        legend.direction = "horizontal",
        title = element_text(size = 13)) +
  ggtitle("A. Coefficient Plot")
fig_sur




#### Simulations

data20 <- read.csv(here("Data", "AME_Bootstrap_CountyLevel_20Shock.RData"))
data40 <- read.csv(here("Data", "AME_Bootstrap_CountyLevel_40Shock.RData"))
data60 <- read.csv(here("Data", "AME_Bootstrap_CountyLevel_60Shock.RData"))
data80 <- read.csv(here("Data", "AME_Bootstrap_CountyLevel_80Shock.RData"))
data100 <- read.csv(here("Data", "AME_Bootstrap_CountyLevel_100Shock.RData"))
data120 <- read.csv(here("Data", "AME_Bootstrap_CountyLevel_120Shock.RData"))


data20 <- data20 %>% 
  mutate(shock = 20)
data40 <- data40 %>% 
  mutate(shock = 40)
data60 <- data60 %>% 
  mutate(shock = 60)
data80 <- data80 %>% 
  mutate(shock = 80)
data100 <- data100 %>% 
  mutate(shock = 100)
data120 <- data120 %>% 
  mutate(shock = 120)

plot_data <- rbind(data20, data40, data60, data80, data100, data120)
glimpse(plot_data)


fig_sim <- plot_data %>% 
  mutate(party = factor(party, levels = c("Republican",
                                          "Democratic",
                                          "Abstention"),
                        labels = c("Republican",
                                   "Democratic",
                                   "Abstention"))) %>% 
  mutate(shock = factor(shock, levels = c(20, 40, 60, 80, 100, 120),
                        labels = c(20, 40, 60, 80, 100, 120))) %>% 
  ggplot(aes(shape = shock)) +
  geom_hline(yintercept = 0, colour = "gray27",  lty = 2, linewidth = 0.75) +
  geom_pointrange(aes(x = party, y = mean, ymin = lb,
                      ymax = ub), colour = "black", fill = "black",
                  lwd = 0.6, position = position_dodge(width = .5)) + 
  scale_shape_manual(name = "Shock in
Death Rates:",
                     values = c(22, 21, 24, 25, 3 , 4)) +
  scale_y_continuous(breaks = c(-0.02, -0.01, 0, 0.01, 0.02, 0.03), lim = c(-0.024, 0.0375)) +
  theme_bw() + #ylim(-0.025, 0.02) +
  labs(y = "Difference in Predicted Share of Voters") +#, 
  #   caption = "Effect of +1 Standard Deviation Change in Unemployment Rate") +
  xlab(NULL) +
  theme(axis.text.x = element_text(size = 13#, 
                                   #hjust = 1, 
                                   #angle = 15
  ),
  axis.text.y = element_text(size = 12),
  axis.title.y = element_text(size = 13),
  legend.text = element_text(size = 13),
  legend.title = element_text(size = 12),
  legend.position = "right",
  legend.direction = "vertical",
  title = element_text(size = 14)) +
  ggtitle("B. Simulated Effect of a Death Rate Shock")
fig_sim



multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  if (numPlots==1) {
    print(plots[[1]])
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

multiplot(fig_sur, fig_sim, cols = 2)






